Hydrodynamics from dynamical 
non-equilibrium MD. 

Sergio Orlandini* ^ Simone Meloni^' **'^ and Giovanni Ciccotti*'^ 

* Dipartimento di Chimica, Universita ''Sapienza'\ Rle A. Mow 5, 00185 Roma, Italy 
^Consorzio Interuniversitario per le Applicazioni di Supercalcolo Per Universita e Ricerca 
(CASPUR), Via dei Tizii 6, 00185 Roma, Italy 
""^School of Physics, Room 302 UCD-EMSC, University College Dublin, Belfield, Dublin 4, Ireland 
^School of Physics, Room 302B UCD-EMSC, University College Dublin, Belfield, Dublin 4, 

Ireland 

^ Dipartimento di Fisica, Universita "La Sapienza'' and CNISM, Piazzale Aldo Moro 5, 00185 

Roma, Italy 

Abstract. 

We review a dynamical approach to non-equilibrium MD (D-NEMD). We show how, using a 
proper simulation setup, is possible to treat interesting cases in which the initial condition is a 
stationary non-equilibrium state produced by a suitable dynamical system. We then extend the class 
of non-equilibrium phenomena that can be studied by atomistic simulations to the case of complex 
initial conditions consisting in assigning a macroscopic value of a scalar or vector observable or a 
field. We illustrate the functioning of this method by applying it to the relaxation of an interface 
between two immiscible liquids. We have shown that our method generate unbiased results while 
this might not be the case for the often used short time average approach. 
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1. INTRODUCTION 

By time averages, the macroscopic properties, that are the object of equilibrium statis- 
tical mechanics, emerge from the microscopic interactions among the elementary con- 
stituents of a macroscopic body in a natural way. This is the Boltzmann point of view put 
in practice by molecular dynamics (MD) simulations. The situation is much less settled 
in the case of non-equilibrium statistical mechanics, except for the linear response to an 
external field [1,2] and for stationary non-equilibrium situations [3, 4, 5, 6] (where the 
use of time averages is still justified). In the first case, the Kubo formula reduces the 
determination of the linear response to the calculation of equilibrium time correlation 
functions, which are easy to sample accurately by MD simulations; in the second case, 
time averages can still be used to replace ensemble averages over the unknown ensemble 
corresponding to the steady state (provide the system remains ergodic). Thus, for both 
cases, with a proper setup the physical properties of the system can be computed. 

A MD method to compute the statistical properties of a non-equilibrium non- 
stationary systems has been proposed quite some time ago [7]. with this method 
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rigorous ensemble averages can be obtained by following the so called dynamical 
approach to non-equilibrium molecular dynamics (D-NEMD) [7, 8, 9]. This method, 
which is just the generalization and the numerical implementation of the Onsager princi- 
ple of regressive fluctuations [10, 11], tells that the time dependent statistical average of 
a microscopic observable is computed by taking an average over the initial ensemble of 
the observable evolved in time under a perturbed dynamics. When the initial ensemble 
corresponds either to a system in equilibrium or to the stationary non-equilibrium state 
of a given perturbation, it can be easily sampled by MD. The Onsager- Kubo relation 
(see below) is thus exploited to perform non-equilibrium MD ensemble averages. 
That involves choosing independent configurations of the system in a steady state, 
performing the perturbed time evolution of every independent initial configuration and 
computing averages over the set of these different trajectories [7, 8, 9]. However, when 
the initial distribution corresponds to the equilibrium distribution of a system submitted 
to a macroscopic constraint, the prescription above must be complemented by a method 
which allows to sample the conditional Probability Density Function (PDF) associated 
with this condition. In this paper we shall illustrate how to deal with both cases: i) initial 
conditions corresponding to a stationary non-equilibrium distribution [12], and ii) initial 
conditions corresponding to the conditional equilibrium distribution of a system subject 
to an external macroscopic constraint, be it of scalar, vectorial or field nature [13]. This 
latter method allows to treat constraint of both scalar or field nature. We will achieve this 
objective by combining D-NEMD with restrained MD. The latter is used to perform the 
conditional average satisfying the macroscopic constraint mentioned above then, like 
in standard D-NEMD, the Onsager-Kubo relation is used to compute time-dependent 
statistical average of the relaxation from this state. It is worth to mention that the scalar 
and vectorial cases have been already implicitly solved by the Blue Moon approach in 
Ref. [14]. However, using the Blue Moon approach for vector or, even worse, field-like 
constraints, especially in molecular systems in which constraints are also used for 
imposing the molecular geometry, might result inconvenient in practice. 

As an illustration of the case of stationary non-equilibrium initial conditions we will 
study the transient hydrodynamical behavior in the formation of convective cells within 
a two-dimensional fluid system subject to orthogonal thermal gradient and gravity field. 
We follow the time-dependence of the density p(x;^), velocity \{x;t) and temperature 
T{x\t) fields induced adding a gravity to a system in stationary conditions under the 
thermal gradient. 

To illustrate our method to sample from a microscopic system under macroscopic 
constraints, we study the hydrodynamic relaxation to the equilibrium of the interface 
between two immiscible liquids. We compute the time evolution of the difference of 
density fields of the species A and B (Ap(x;^) = p^{x\t) — p^{x\t)) and the associated 
velocity field (v^ (x; t)) and show the usefulness of our approach that does not rely on the 
separation of timescale of atomistic and hydrodynamical processes, as it is the case for 
the often used method of local time average. 

The paper is organized as follows. In Sec. 2 we describe the theoretical background 
of D-NEMD and our method for sampling complex initial conditions. In Sec. 3 we 
present an application of D-NEMD to the study of the formation of convective cells in a 
system subject to thermal gradient and gravity field. In Sec. 4 we discuss an application 
of our combined restrained MD + D-NEMD method. Finally, in Sec. 5 we draw some 
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conclusions. 



2. THEORETICAL BACKGROUND 

In D-NEMD we consider a classical system with N particles in d spatial dimensions. 
Let r/ and p^- be the position and the momentum of the i-th particle, respectively, and 
r = {r/, p^} be a point in phase space. The Hamiltonian governing the system is H(r, t), 
which we assume to be explicitly time dependent (the generalization to non-Hamiltonian 
systems could be easily worked out [16, 17, 18]). 

In statistical mechanics an observable, including a macroscopic field (see below for 
its definition), is obtained as an ensemble average in phase space of the corresponding 
microscopic observable: 

o{t) = J drd{r)w{r,t) = (d(r)w(r,o> (i) 

where 0{t) is the macroscopic observable and 0(r) is the corresponding microscopic 
observable. w(r^t) is the (time dependent) phase-space PDF. As a consequence of the 
conservation of probability, w{r,t) obeys the Liouville equation: 

^ = -VrlMm (2) 

{w{r,t),H{r,t)} = -iL{t)w{r,t) 

where {•, •} is the Poisson bracket and L{t) is the Liouville operator. A formal solution 
of Eq.(2) is: 



w(r,o = sHt)w{r,o) o) 

where S'^{t) is the adjoint time evolution operator of the dynamical system. 

We now consider the time-evolution of a non-time-dependent phase-space observable: 



dt 

which has the following formal solution 



"^^^^^ = {o(r),i/(r,0} = /L(Oo(r) (4) 



o(r(0) = s{t)6{r{o)) (5) 

where S{t) is the time-evolution operator. 

By combining Eq.(l) with Eq. (3) and (5) we get the Onsager-Kubo relation: 



Hydrodynamics from dynamical non-equilibrium MD. 



January 9, 2013 



3 



o{t) = <d{r)s\t)w{r,o)> (6) 
<5(0o(r)w(r,o) > 

The meaning of the above equation is that the ensemble average of the microscopic 
observable 0(r) over the time-dependent PDF w(r,t) at time t is the same as the 
ensemble average of the microscopic observable at the point r(^), corresponding to the 
evolution in time of the initial phase-space point r(0), averaged over the PDF at time 

^ = o(w(r,o)). 

In the standard D-NEMD approach we assume that w(r, 0) corresponds to a stationary 
condition, then can be sampled via an MD simulation possibly governed by the Hamil- 
tonian Ho{r) (more general stationary conditions can be constructed, see below and in 
Ref. [12]). Then we can evolve initial configurations taken from that trajectory with the 
dynamics generated by H(r,t). Along these paths we compute the microscopic observ- 
able d(r{t)). The time-dependent behavior of the macroscopic observable 0{t) is the 
ensemble average of d(r{t)) over all the trajectories originated from each of the initial 
states (see Eq. (6)). An example of how to generate stationary non-equilibrium condi- 
tions might be by putting the sample in contact with two thermal reservoirs at different 
temperature at the border of the simulation box. Using this setup it is possible by MD to 
sample an initial stationary PDF corresponding to a system under thermal gradient. 

We will now generalize the standard D-NEMD approach to the case of an initial PDF 
associated with a macroscopic constraint, i.e. the conditional PDF to find the system at 
the point F in phase space given that a microscopic field ^(x|r), function of the atomic 
position only, at the time ^ = is equal to F* (x) (F(x|r) = F* (x)). Our method a fortiori 
applies to the case in which the observable F is a scalar. For sake of clarity we first 
introduce our method assuming that F is scalar and then explain how to extend it to the 
case of space-dependent fields. Let us start recalling that the conditional PDF mentioned 
above reads: 



.(r,F = F.) = !«r)^ (7) 

where F(r) is the microscopic observable associated with the macroscopic observable 
F, ^ is the partition function associated with w(r) and Pf(F*) = / drw(r)5{F(r) — 
F*) is the PDF to find the system in the state F ^F"". 

We propose to sample this conditional PDF by the biased MD simulation governed 
by the following Hamiltonian: 

H{T)=Ho{T) + ^-{F{v)-F*f (8) 
where A: is a tunable parameter. The PDF sampled by this biased MD is: 
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exp[-i6(//o(r) + |(F(r)-F*)2)] 
/jrexp[-j8(H(r,0) + |(F(r)-F*)2)] 



By recalling that exp(-(j - jla^) / \flno ju), in the limit of ^ ^ oo, we 

see that Eq. (9) goes to 

w rn ^ exp[-i6(//o(r)]g(F(r)-F-) 
^ /jrexp[-^(/fo(r)]5(F(r)-F*) 

Multiplying and dividing by ^ it is apparent that Wj^(r) w(r|F = F*). 

We now extend our discussion to the case in which the condition (and the observ- 
ables of interest) is with respect to a macroscopic field, rather than to a scalar. Generally 
speaking, in statistical mechanics the microscopic observable associated to a macro- 
scopic field is defined as [21]: 

F(x,r) = J^^Kr)5(x-rO (11) 

/=i 

where ^/(F) is the microscopic property under consideration relative to the particle 
/ and X is a point in the ordinary space. Eq. (11) means that only atoms at the 
point r/ = X contribute to the field at that point. Examples of fields of interest for the 
problem presented in Sec. 3 and 4 are the fields density (p(x,^)), velocity (v(x,^)) and 
temperature (r(x,^)). These fields at time t are given by: 

p(x,0 =< J2M/5(x-rXr,0 > (12) 

<if^.M(;;-;.Mr.o> 

T(^t) = — <I^li^i^-^i)[Pi-^i'^MM^,t)> ^^^^ 

3kB P(x,0 

where fii is the mass of the i-th atom. 

Eq. (7) can be easily extended to the case in which the condition is a field: 

Now w[r|F(x)] is functional of the field F{x) and the notation 5(F(x|r) — F*(x)) = 
rixGSt^ 5(F(x,r) — F*(x)) indicates that the delta function is considered acting all over 
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the X space. Of course, this definition of the conditional probabihty density functional 
is of no practical use in simulations. We overcome this problem by introducing a 
discretization {xa}a=i,m 

of the 91^ 

space, where m is the number of points over which 
the space is discretized. Consistently with this discretization the microscopic observable 
fields is defined as the average over the cells around around each point Xo^: 

F{xa,r) = ^[ Jx£^Kr)5(x-r,-) = ^£^/(r) (i6) 

where the sum in the r.h.s. runs over the atoms belonging to the cell around the point x^. 
However, this definition is not suitable for our restrained MD method as it might produce 
impulsive forces on the atoms moving from one cell to another (see below). Therefore 
we smooth the central term of Eq. (16) by replacing 5(x — r/) with the (normalized) 
gaussian g(x;r/,cj) = exp[— (x — r/)^/2(7^]/v/27rcT. With this approximation, Eq. (16) 
becomes: 



^(xa,r) = 



(17) 



where erf{c^G) = fl^dx g(x;0, c) is the error function, and Xa^x ^a,x ^ 
component of the upper and lower limit of the orthorombic cell around the point Xa. 
The conditional PDF associated to the field F{x) on the discrete representation of the 
space is: 



w(r|{F(x«,r) -F (x«)}«=,.) (18) 

where /V({^*(xa)}a=i,m) is the joint probability that F(xi,r) = /^*(xi) 
F(x^,r) = F*(x^). This conditional PDF can be sampled by a MD governed by 
the following Hamiltonian 



^ k 

H{r)=Ho{r)+ ^ -{F{xa,r)-F*ixa)f (19) 
a=l ^ 

which is the straightforward extension of the Hamiltonian of Eq. (19) to the case of a 
macroscopic field in the discrete space approximation. Then, as in standard D-NEMD 
approach, we evolve a set of initial configurations taken from the trajectory above with 
the dynamics generated by H(T^t). Along these paths we can compute the microscopic 
observables and calculate the ensemble average over all the trajectories originated from 
each initial state (Eq. (6)). 

Before closing this section it is worth to mention that the problem of sampling 
initial conditions consistent with a macroscopic constraint was already implicitly solved 
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by the Blue Moon ensemble (see Refs. [14] and [15]). Blue Moon does not sample 
directly the conditional PDF w(r|F = F*) (or w (r|{F(xa,r) = F%Xa)}a=\,m) for the 
vectorial case) but rather it samples the constrained PDF in configurational space wp* (r) . 
The relation between the Blue Moon ensemble average and conditional average in the 
configurational space is given by (see Ref. [15]): 



<d(rn/vmm\>B. 

< l/\/|det[C(r)]| >BM 

where 0(r) is a microscopic observable andC/j (r) = Vcj/(r)M~^Vcjy(r), being cj/(r) the 
(vectorial) condition F(r) = F* or any other constraint imposed on the system, in par- 
ticular constraints imposed for modeling (partly) rigid molecules [19], and M the mass 
matrix (Mij = liidif). The F — F"" and BM indexes denote that the ensemble averages are 
taken according to the conditional or constrained (Blue Moon) PDF, respectively. If the 
observable of interest depends on the phase space (0(r)) rather than on the configura- 
tional space it is possible to extend the validity of Eq. 20 by generating the momentum 
component of the PDF from a suitable Maxwellian distribution. 

Two comments are in order concerning the restraint method for sampling the complex 
initial condition described in this paper versus the approach based on the Blue Moon 
ensemble. First of all, depending on the macroscopic condition and the molecular con- 
straints, calculating the unbiasing term 1 / ^|det[C(r)]|/ < 1/ -^|det[C(r)]| >bm might 
be complex. Moreover, while the restraint approach can be combined with Monte Carlo 
simulation when the microscopic observable connected to the macroscopic condition is 
not analytical (see Ref. [20] for a detailed description of this approach), the same cannot 
be done with the Blue Moon approach. 

In the following two sections we show two applications of the D-NEMD described 
above. We first present the case of an initial condition corresponding to a stationary non- 
equilibrium system. In particular we study the transient state of formation of convection 
cells when a gravity force is added to a system subject to a thermal gradient (Sec. 3). 
Then we illustrate an application of the restrained method to sample initial conditions 
corresponding to a system subject to macroscopic constraints by studying the relaxation 
to the equilibrium of the interface between two immiscible liquids (Sec. 4). 



3. FORMATION OF CONVECTIVE CELLS 

3.1. Computational Model and Setup. 

A fluid system consisting of N = 5041 particles is contained in a two- dimensional 
box in the {xz} plane, with the gravity force directed along the negative verse of the z 
axis. The particles interact via a WCA (Weeks-Chandler- Andersen) potential 

r 4e[(c7/r)i2-(c7/r)6],Vr<2i/6a 



Hydrodynamics from dynamical non-equilibrium MD. 



January 9, 2013 



7 



where £ and a are the usual Lennard- Jones parameters. This potential is a purely repul- 
sive potential obtained from the Lennard- Jones 12-6 potential, truncated in its minimum 
(so that the force is continuous), and shifted (so that the potential is continuous). Each 
particle has mass /i. Hereafter, we will use reduced units putting the typical scale of 
energy, mass, and length equal to unity, i.e., 8 = 1, jU = 1, and CJ = 1. Times are then 
measured in units of ^ jia^ /e. 

We want to study the system in physical conditions that allow the formation of a 
convective cell, i.e., when gravity and a thermal gradient are present. Moreover, we want 
to analyze the transient evolution to the formation of the steady-state roll. For this we 
can take as initial condition of the system the steady state under the effect of a thermal 
gradient and then study the dynamical response of the system to the ignition of gravity. 
As far as the physical setting is concerned, we take the thermal gradient orthogonal to 
the gravity force. This allows a straightforward application of the D-NEMD technique. 
At the top and at the bottom of the box there are repulsive walls to avoid particles from 
drifting downwards under the effect of the gravity force. The thermal reservoirs which 
produce the thermal gradient are realized as two stripes such that the components of the 
velocity of each particle located in one of these stripes are sampled from a Maxwellian 
distribution at the temperature of the wall. We assume periodic boundary conditions at 
the thermal walls located at the two lateral sides of the box, i.e., a particle can move from 
the hotter to the colder reservoir. To avoid that particles near a thermal reservoir interact 
simultaneously with both reservoirs, we chose the thickness of the reservoir — 1.68, 
larger than the cutoff of the WCA interaction. When the system is in equilibrium, each 
reservoir contains roughly 100 particles on average. We chose the temperature of the 
colder reservoir as Ti = 1.5 and a theoretical thermal gradient |Vr| =0.1, so that the 
hotter reservoir has a temperature T2 — |Vr|L + Ti — 9.9. The gravity force used is 
g = 0.1 in Lennard- Jones units. In SI units, taking Ar as a reference fluid, we have 
g = 7 X 10^2 m/s^ T2 = 1 196 K, and Ti = 179.7 K. The box length is L = 2.89 x 10"^ m, 
so that the thermal gradient is \VT\ = {T2 -Ti)/L = 3.52 x 10^^ K/m. External fields of 
this strength are necessary to reach a sufficient signal to noise ratio in a small system. The 
density, velocity and temperature fields are computed on a 15 x 15 discretization of the 
space. Finally, The averages for these simulations are performed over an ensemble 
made of 1000 copies of the system 



3.2. Results and discussion. 

We follow the evolution of the system after the igniton of the gravity field by mon- 
itoring the time evolution of the temperature and density fields in some characteristic 
cells (see Fig. 1). These fields become stationary at ^ = 250. During the transient, both 
temperature and density oscillate with nearly the same period T = 18 in all cells. The 
phase of these oscillation is constant in all cells at the bottom of the box, and it is the 
same for p and T. The same feature holds for the phase in the cells at the top of the box, 
which is, however, opposite to the phase in the cells at the bottom. During the transient, 
an increase (decrease) in T corresponds to an increase (decrease) in p in the same cell, 
while, in the stationary state at large time, a temperature lower (higher) than the initial 
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value is associated with a density higher (lower) than the initial value, as expected. 
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FIGURE 1. Density (left) and temperature (right) fields as a function of time in selected cells near the 
hot reservoir (top panel), in the middle (middle panel) and near the cold reservoir (bottom panel). For a 
given distance from the reservoirs, three cells at different height are shown. 

This behavior can be followed by plotting the transient velocity field at every quarter 
period of the density and temperature oscillations (Fig. 2). At ^ = t/4 (Fig. 2/a) the 
velocity field points downwards as a consequence of the ignition of gravity at f = 0. 
Kit — xjl (Fig. 2/b), the velocity field is almost null, and the fluid is almost at rest 
on average. At f = 3t/4 (Fig. 2/c) the fluid is expanding against the gravity force 
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and the velocity field is directed upward, as a reaction to compression. Once again, 
in correspondence to the subsequent relative maxima or minima (Fig. 2/d) the fluid 
is nearly at rest. During the following compressions and expansions the flux becomes 
localized, respectively, near the cold the hot reservoirs (see Fig. 3). In correspondence 
of the next relative maxima and minima the fluid is no longer at rest, and instead it begins 
to support a convective flow. The sequence of compressions and expansions strengthens 
the flux, which becomes stable at ^ = 250. The final shape of the convective pattern is 
symmetric (see Fig. 4). 
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FIGURE 2. Local velocity field averaged over 1000 independent initial configurations during the first 
oscillation of temperature and density fields, i.e., at (a) ^ = 4.5 ~ t/4, (b) ^ = 9 ~ t/2, (c) ^ = 13 ~ 3t/4, 
and (d) ^ = 18 ~ T. Hot reservoir on the left side of the box, cold reservoir on the right side of the box. 
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FIGURE 3. Local velocity field averaged over 1000 independent initial configurations during the 
first oscillation of temperature and density fields, i.e., at (a) ^ = 22.5 ^ 5t/4, (b) r = 27 ~ 3t/2, (c) 
r = 31.5 ~ It/4, and (d) t = 36 ^ 2t. Same conditions as in Fig. 2 



4. RELAXATION OF THE INTERFACE BETWEEN TWO 

IMMISCIBLE LIQUIDS 



4.1. Computational Model and Setup. 



We now illustrate the method described in Sec. 2 for sampling complex initial con- 
ditions by studying the relaxation to equilibrium of the interface between two immis- 
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FIGURE 4. Local velocity field averaged over 1000 independent initial configurations at f = 250, i. e. 
at the end of simulation. Same conditions as in Fig. 2 



cible liquids. We define the interface between the two Hquids as the surface at which 
Ap(x) = p'^(x) — p^(x) = 0, where p'^(x) and p^(x) are the densities of the Uquids A 
and B (see Eq. (12)), respectively. We start from the interface defined below and follow 
the isosurface S{t) = {x : Ap{x,t) = 0} of the fields Ap{x,t) and V^(x,f) till equilib- 
rium. The initial conditional PDF is sampled using the method described in Sec. 2 with 
the restraint that Ap(xa) = in the cells, centered around the points x^, through which 
passes the following surface: 

where ^ is the amphtude of the curved surface and {L^}^=i^i, is the length of the 
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simulation box along the ;f-th cartesian direction. The terms and — ^ are added to 
place the interface at the centre of the simulation box. The condition above is the discrete 
counterpart of the continuos condition Ap(x) = 0, V x G 5. We do not impose any other 
condition on the density. However, we prepare the system such that all the particles on 
one side of the interface are of one kind, say A, and of the other kind on the other side, 
say B. Since we apply periodic boundary conditions along all the cartesian directions, 
we have a second flat interface at the beginning/end of the simulation box. 

The sample used in our simulation consists of 171,500 particles: 88,889 of type 
A and 82,611 of type B. Two particles of the same type interact via Lennard- Jones 
potential i/^{r) = u^^{r) — Ae [(cj/r)^^ — (cj/r)^] , while two particles of different type 
interact via the repulsive potential u^^{r) = Ae [(a/r)^^] , where £ and a are the usual 
Lennard- Jones parameters. The simulation box is a parallelepiped of size ^ 45 x 45 x 90 
in reduced units (average density p = l.OlAparticles / G^). In the restrained MD the 
temperature of the sample is kept fixed at l.Ss/kB (ks Boltzmann constant). This density 
and temperature are in the fluid domain of a pure Lennard- Jones system. The ordinary 
space is discretized in 5488 points (a 14 x 14 x 28 grid) and each cell contains, on 
average, ^ 30 particles. 

The system is prepared by thermalizing a sample of pure type A Lennard- Jones 
particles at the target temperature and density, and then transforming those particles 
belonging to the cells on one side of the interface in type B particles (see Fig. (5)). In 
the cells belonging to the interface 5(x) only half of the particles are transformed from 
A to B, so as to have Ap(x) = in these cells. The system is then thermalized with 
the restraint on the Ap(x) for 1.6 x 10^ timesteps. Such a very long run is needed to 
relax the gradient of temperature formed when the nature of particles on one side of 
the interface is changed from A to B (immediately after the A-to-B transformation the 
interfaces - the curved one and the flat one due to the periodic boundary conditions - 
due to the strong repulsive forces among particles of different type, are warmer than the 
bulk). The timestep used in this and the next phase (restrained MD runs) is 4.56 x 10""^ 
LJ time units, which is one order of magnitude smaller than the typical timestep for 
simulation of Lennard- Jones systems. This very short timestep is required by the stiff 
force associated to the restraint. After this relaxation, a 10^ timestep long restrained MD 
is performed along which, at regular intervals of 25,000 timesteps, we collect 40 initial 
positions and velocities for the second step of the D-NEMD procedure. The atomic 
configuration corresponding to one snapshot of this trajectory is shown in Fig. (5). 

In Fig. (6) is reported the Ap(x,^) field on the points {xa}a=i,m together with the 
isosurface Ap (x) = obtained as a linear interpolation of the value of the field on the 
grid points. This figure shows that the interface is rather sharp, involving typically one 
or two shells along the direction orthogonal to it. 

From each of these initial conditions we start 25,000 timestep long unrestrained 
MD simulations along which we compute the microscopic observables of interest. By 
averaging over the (40) initial conditions we get p^{x^t) and v'^(x,^). 



Hydrodynamics from dynamical non-equilibrium MD. 



January 9, 2013 



13 




FIGURE 5. One atomic configuration extracted from the restrained MD used for sampling the initial 
conditional PDF. In blue particles of type A, in red particles of type B. 



4.2. Results and discussion. 

In this section we present our results obtained with the restrained MD method and 
compare them with those obtained by computing the relevant fields along one unre- 
strained trajectory started from a configuration sampled from the restrained dynamics. 
This second type of simulation, often combined with a "local time average" (i.e. aver- 
aging over a small time- window centered at the current time), are used to study hydro- 
dynamical phenomena by atomistic simulation. We show that the fields obtained from 
this latter type of simulation violate some of the properties of the hydrodynamical fields 
associated to the process under investigation while our D-NEMD approach does not. 

Let us start by analyzing the surface S{t) = {x : Ap(x,f) = 0}. In Fig. (7) it is shown 
a series of snapshots of the S{t) surface. First of all we remark that all along its evolution 
the surface satisfies the symmetry of the problem, i.e. it is symmetric with respect the 
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FIGURE 6. Ap(x) field over the grid points {xa}a=i,m- The color of each point depends on the Ap(x) 
on that point. Intense red means that all the particles in the cell are of type A and intense blue means that 
they are all of type B. Intermediate colors indicate that the cell contains both types of particles (white 
correspond to 50% of each type of particle). The curved interface represent the isosurface Ap(x) = as 
obtained from the linear interpolation of the density on the grid points. The second interface due to the 
periodic boundary conditions on the long edge of the simulation box is not shown. 



{yz} reflection plane passing by the middle of the simulation box and it is translationally 
invariant along the y direction. The tiny bumps on the S{t) surface are due to the limited 
number of initial conditions (40) used for computing Ap(x,f). In the limit of an infinite 
number of such initial conditions the surface would be completely smooth, as expected 
and predicted by classical hydrodynamics. The relaxation from the initial curved surface 
to the final flat surface takes approximately 20,000 timesteps which, if we chose for the 
Lennard-Jones parameters the values a = 3.405 A and £ = 0.01032 eV (suitable for 
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modeling Ar) and an amplitude of the initial interface ^ = 50 A, corresponds to a 
maximum value of the field v(x) of ~ 80 m/^". 




FIGURE 7. Snapshots of the interface S{t). The arrows on the first snapshot show the direction of 
evolution of the interface: the center and the extreme of the interface move in opposite directions. The 
field Ap(x,^) is also shown on one {xz} plane by adopting the same colorcoding of Fig. 7 

We now move to the analysis of the velocity field, focusing on the velocity field of 
only the chemical species A: 



v^(x,0 = (23) 



Irolt'iP/(ro)n3^ 


1 


erf{xa,x - r,- ;t(ro), ct) - erfix^ ^^ - r,- ;^(ro), cr) 




erf{xa,x - ;^(ro), ct) - erf^x^ ^ - r,- ;^(ro), o) 





where the sum runs only over the atoms of type A, and it is implicitly assumed 
that and r/ are taken at the time t starting from the initial condition Fq (see Eqs. (17) 
and (13)). The sum Y.Tq ^^^^ the initial conditions along the restrained MD. We 
consider the field v"^ (x, ^ ) as, due to the conservation of the total momentum and the fact 
that the initial total momentum was set to zero, the total field, i.e. those including A and 
B specie, is, on average, zero. In the left column of Fig. 8 is shown the velocity field 
v'^(x,^) at various times. This field is computed only on the grid points corresponding to 
cells that contains at least one particle of type A. This fact makes the field "noisy" (large 
values of the field rapidly changing orientation) close to the A/B interface, where the 
cells contain less A particles, and therefore the average of the atomic velocities over the 
particles in the cell (see Eqs. (16) and (17)) is less effective in smoothing the field. This 
effect is reduced by making larger the number of initial configurations used to perform 
the ensemble average over the initial conditional PDF. As a first remark, it is worth to 
mention that relatively few cell layers nearby the interface are involved in the relaxation 
process. In fact, already 5 — 10 cells far from the interface the velocity field is essentially 
zero at any time during the relaxation. Coming to the hydrodynamical process producing 
the relaxation of the interface, we notice that initially (see panel 1 of Fig. 8) the velocity 
field at the top of the interface is pointing downward while at bottom it is pointing 
upward. After some time this field stabilizes into a double symmetric roll, one rotating 
clockwise and the other one rotating counter clockwise, both starting at the top of the 
interface and ending at its bottom (see panel 2 of Fig. 8). Overall, this velocity field 
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produce the phenomenon of pushing up the side of the interface and puUing down the 
center as shown in Fig. 7. The relaxation of the interface follows this mechanism almost 
till the end of the process. In fact, in panel 3 of Fig. 8 we see that still after 45.6 LJ time 
units (10^ timesteps), when the interface is almost flat, the double roll is still present. 
Eventually, after 114 units of time (2.5 x 10^ timesteps) the interface is completely flat 
and the field is null everywhere (panel 4). 



v^(x,i) v"(x,i) 



A 



■■111 3^ 




1;-"^ t= 0.456 

CJI^ / \ 



TIT 



1) . 



I / 



t= 4.56 



t= 45.6 



' f t \ 



t= 114 



FIGURE 8. Snapshots of the \^{x,t) (left, see Eq. (13)) and \^{x,t) (right, see Eq. (24)) fields on the 
grid points belonging to one {xz} plane at various times. 

It is very interesting to compare the V^(x,f) field as obtained from the D-NEMD 
simulation with the instantaneous field v^(x,^) defined as 









-^rf{x^^^-ri^X^o)\^ 








-erfi^a,x-^hX^^)\ 



(24) 



Few snapshots of this field are shown in Fig. 8. First of all we notice that the interface 
relaxation process occurs via the formation of a clockwise roll, which is initially at the 
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top of the interface (panel 1 of Fig. 8) and then move toward the bulk (panel 2 of the 
same figure). This roll is stable all over the duration of the simulation (see panel 3) and 
eventually disappear when the equilibrium is reached (panel 4). The shape of the field 
v'^(x,^) contrasts with the symmetry imposed on the problem, which implies a {yz} 
mirror plane passing through the middle of the simulation box. This problem cannot be 
solved by the "local time average" that is often used in simulation of hydrodynamical 
processes by atomistic simulations. In fact, as mentioned above, the clockwise vortex is 
very stable and a local time average will not restore the proper symmetry expected for 
this field. This problem illustrates that a proper statistical average is needed in order to 
compute by atomistic simulation hydrodynamic fields as otherwise a unlucky choice of 
the initial conditions can produce unphysical results. 

5. CONCLUSIONS. 

In this paper we have reviewed the dynamical approach to non-equilibrium MD. We 
have shown that using a proper simulation setup it is possible to treat interesting cases 
in which the initial condition is either a stationary non-equilibrium condition or a 
constrained equilibrium consistent with the value of a macroscopic scalar or field-like 
observable. We illustrated the functioning of the method by applying it to two cases: 
the establishing of convective cells and the relaxation of an interface between two 
immiscible liquids. We have shown that our method generates rigorous time-dependent 
non-equilibrium averages, while the method of local time average, often used to simulate 
hydrodynamical processes from atomistic simulation, can, sometimes, fail. 

Our conclusion is that the method is ready for challenging applications. Work is in 
progress in this direction. 

ACKNOWLEDGMENTS 

S. M. and G. C. acknowledge SFI Grant 08-IN.1-I1869 for the financial support. S.O. 
acknowledges SimBioMa for the financial support. Finally, the authors wish to acknowl- 
edge the SFI/HEA Irish Centre for High-End Computing (ICHEC) for the provision of 
computational facilities. 

REFERENCES 

1. B. J. Alder, and T. E. Wainwright, Phys. Rev. A 1, 18 (1970). 

2. T. Wainwright, B. Alder, and D. Gass, Phys. Rev. A 4, 233 (1971). 

3. A. W. Less, and S. F. Edwards, /. Phys. C 5, 1921 (1972). 

4. E. Gosling, I. McDonald, and K. Singer, Mol. Phys. 26, 1475 (1973). 

5. W. Ashurst, and W. Hoover, Phys. Rev Lett. 31, 206 (1973). 

6. W. G. Hoover and W.T. Ashurst, Adv. Theor. Chem. 1, 1 (1973). 

7. G. Ciccotti, and G. Jacucci, Phys. Rev. Lett. 35, 789 (1975). 

8. G. Ciccotti, G. Jacucci, and I. R. McDonald, /. Stat. Phys. 21, 1 (1979). 

9. G. Ciccotti, C. Pierleoni, and J. P. Ryckaert, in Microscopic Simulations of Complex Hydrodynamic 
Phenomena, Plenum, New York, 1992. 

10. L. Onsager, Phys. Rev. 37, 405 (1931). 



Hydrodynamics from dynamical non-equilibrium MD. 



January 9, 2013 



18 



11. R. Kubo, /. Phys. Soc. Jpn. 12, 570 (1957). 

12. M. L. Mugnai, S. Caprara, G. Ciccotti, C. Pierleoni, and M. Mareschal, /. Chem. Phys. 131, 064106 
(2009) 

13. S. Orlandini, S. Meloni, and G. Ciccotti, in preparation 

14. E. A. Carter, G. Ciccotti, J. T. Hynes and R. Kapral, Chem. Phys. Lett. 156, 472 (1989) 

15. G. Ciccotti, R. Kapral, and E. Vanden-Eijnden, Chem. Phys. Chem. 6, 1809 (2005) 

16. M. E. Tuckerman, C. J. Mundy and G. J. Martyna Europhys. Lett. 45, 149 (1999) 

17. M. E. Tuckerman, Y. Liu, G. Ciccotti, and G. J. Martyna, /. Chem. Phys. 115, 1678 (2001) 

18. C. Hartmann, C. Sch§tte and G. Ciccotti, J. Chem. Phys. 132, 1 1 1 103 (2010) 

19. G. Ciccotti and J. R Ryckaert, Comp. Phys. Rep. 4, 347 (1986) 

20. G. Ciccotti and S. Meloni, submitted to Phys. Chem. Chem. Phys. 

21. J. H. Irving, and J. G. Kirkwood, /. Chem. Phys. 18, 817 (1950). 



Hydrodynamics from dynamical non-equilibrium MD. 



January 9, 2013 



19 



